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Evidence for chaotic behaviour in pulsar spin-down rates 
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ABSTRACT 

We present evidence for chaotic dynamics within the spin-down rates of 17 pulsars 
originally presented by Lyne et al. Using techniques that allow us to re-sample the 
original measurements without losing structural information, we have searched for ev- 
idence of a strange attractor in the time series of frequency derivatives for each of 
the 17 pulsars. We demonstrate the effectiveness of our methods by applying them 
to a component of the Lorenz and Rossler attractors that were sampled with similar 
cadence to the pulsar time series. Our measurements of correlation dimension and Lya- 
punov exponent show that the underlying behaviour appears to be driven by a strange 
attractor with approximately three governing non-linear differential equations. This 
is particularly apparent in the case of PSR B1828— 11 where a correlation dimension 
of 2.06 ± 0.03 and a Lyapunov exponent of (4.0 ± 0.3) x 10~ 4 inverse days were mea- 
sured. These results provide an additional diagnostic for testing future models of this 
behaviour. 

Key words: chaos — methods: data analysis — stars: kinematics and dynamics — 
stars: rotation — pulsars: general — pulsars: individual: B1828— 11 



1 INTRODUCTION 



Pulsars are spinning neutron stars whose emission is 
thought to be driven by their magnetic fields (see, e.g., 
IManchester fc Taylor! Il977h . Their dynamics exhibit a wide 
ranging degree of stability, with the fastest spinning and 
oldest ' millisecond pulsa r s' gen erally being the most pre- 
dictable. IPetit fc Tavellal (| 19961 ) showed that these millisec- 
ond pulsars can be as stable as atomic clocks on large time- 
scales. 

Because pulsars are so stable, it is surprising when 
we see them misbehave. Departure from normal behaviour, 
characterized by steady emission and rotation, can occur in 
a variety of ways. Some pulsars have nulling events, where 
the emission s eems to turn off for a while and then suddenly 
turn back on <|Backerlll97cl ). Even more extreme behaviour 
can be seen in the intermittent pulsar B1931+24, which be- 
haves like a normal pulsar fo r five to ten days, t hen is un- 
detectable for about 25 days |Kramer et al.ll2006t ). Recently 
even longer-term intermittency (spannin g hundreds of days) 
has been reported in two other pulsars ( Camilo et al.ll2012l . 
Lorimer et. al. 2012 in prep.). Finally, another related class 
of pulsars known as Rotating Radio Transients (RRATs) 
seem to sporadically turn their emission on and off on a 
wide range of time-scales. 

It is even more shocking when we see changes in the 
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dynamics of a pulsar. Though pulsars are expected to grad- 
ually spin-dow n over time , due to an energy loss from mag- 
netic braking (|Goldlll969l ). other changes in the dynamics 
are highly unexpected. These fluctuations often give us in- 
sight to the interior and the environment of a pulsar. One 
such phenomenon, known as a 'glitch', is a sudden discrete 
increase in rotation that has been attrib uted to superfluid 
vorti ces within the interior of the pulsar (|Anderson fc Itohl 

Hazi). 

Since the dynamical fluctuations in pulsars were unex- 
pected, a bulk of the irregulariti es were considered to be 
'timing noise' llHobbs et al.l l2006). When these irregulari- 
ties were observed over a 20 yr time span, larg e time-scale 
fluctu ations in the spin-down rate became clear. ILvne et al.l 
l|2010l ) describe the irregularities as 'quasi-periodic' and were 
able to relate them to the pulse shape. There have been sev- 
eral differe nt processes proposed for these fluct uations from 
precession (| Jones! I20T2} ) to non-radial modes (| Rosen et al.l 
1201 if ) . Yet, the mechanisms that govern these fluctuations 
and their connection to the pulse shape is still a mystery. 
Quasi-periodicities are often a sign of a non -linear chaotic 
system. Previous chaotic studies o n pulsars dHarding et al' 
| l990l ; iDelanev fc Weatheralll 1 19981 : iDeLanev fc Weatherali 
Il999l ) focused on emission abnormalities and timing noise 
for particular pulsars. 

In this paper we search for chaotic behav iour within 
the sp in-down rate of 17 pulsars presented by ILvne et al] 
|20ld ). In Section 2 we give a brief introduction to chaotic 
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(a) 



(b) 



Figure 1. (a) The x component of the Lorenz equations whose 
initial conditions differed by 10 — 3 in x. (b) The difference be- 
tween the two scenarios in (a) over time. An overall increasing 
exponential trend can be seen. 



systems and their behaviours. In Section 3 we present tech- 
niques to form an evenly sampled time series with the same 
structural information as an unevenly measured series. In 
Section 4 we demonstrate methods to search for different 
chaotic characteristics within a time series, and test our al- 
gorithm with known chaotic systems. In Section 5 we discuss 
our results and their implications. The techniques presented 
here are very general and can be used on almost any time 
series. Therefore, we have written this paper explicitly in 
the hope that these techniques will be more approachable, 
and that they will be used more frequently within the pulsar 
community. 



2 CHAOTIC BEHAVIOUR 

In everyday conversation, the word chaos is often inter- 
changeable with randomness, but in dynamical studies these 
are two distinct ideas. Chaos is continuous and determinis- 
tic with underlying governing equations, while randomness 
is more complex and uncorrelated; values at an earlier time 
have no effect on the values at a later time. 

One of the chara cteristics of ch aos has been colourfully 
described by Edward lLoren j l| 19931 ) as the 'butterfly effect'. 
Lorenz asks, 'Does the flap of a butterfly's wings in Brazil set 
off a tornado in Texas?' This is used to illustrate an insta- 
bility, where a system is highly sensitive to initial conditions. 
The consequence is that if there is a small displacement in 
the initial conditions, the difference between the two scenar- 
ios will grow exponentially to cause significant changes at a 
later time. This instability does not arise in linear d ynamics 
and i s a chaotic phenomenon in non-linear systems (jScargld 
ll992T l. We will utilize this chaotic trait in Section \4. 31 

An example of this behaviour is shown in Fig. [1] for the 
Lorenz system of equations: 

x — a(y — x) 

y = x(p-z)-y (1) 
z = xy — f3z. 

Here a, p, and /3 are positive parameters. lLoren 1 (|1963T I 
derived this system from a simplified model of convection 
rolls in the atmosphere. In the original derivation a is the 
Prandtl number, p the Rayleigh number and f3 has no proper 
name but relates to the height of the fluid layer. As for 
the governing variables, x is proportional to the intensity of 
the convective motion, y is proportional to the temperature 
difference of the acceding and descending currents, and z 



is proportional to the distortion of the vertical temperature 
profile from linearity l|Lorenzlll96St ). Since then, the Lorenz 
equations have appeared in a wide range of physical systems. 

It is important to note that there are no analytical so- 
lutions to most non-linear equations. Often to produce a 
solution, as used in Fig. [2j the system is marched forward 
with small enough time steps that enable linear relationships 
to be used to simulate a function. We can see from the gov- 
erning equations that any function that is produced will be 
highly dependent on the other variables. When these time 
functions are plotted with respect to each other, as seen in 
Fig. [2ja), they trace out a rather odd surface. 

The Lo renz equations are not the only system in which 
this occurs. iRosslerl l| 19761 ) was in search of a simpler set of 
equations with similar chaotic behaviour to the Lorenz at- 
tractor. He came up with with the following three equations 
with only one non-linear term zx: 



x = —y — x 
y = x + ay 
z — b — zix — c). 



(2) 



Here a, b, and c are parameters. Again, plotting the variables 
against each other produces a different bizarre surface, seen 
in Fig. 



2(b) 



Though the Rossler equations started solely 
as a mathematical constru ct, analogous behav iour has been 
seen in chemical reactions (lArgoul et al.lll987l ). 

The complex shapes that the dynamics form are known 
as 'strange attractors'. They are called 'attractors' because, 
regardless of the initial conditions, the functions will con- 
verge to a path along these surfaces. They are 'strange' be- 
cause they are fractal in nature. The dimension of the at- 
tractor will be a non-integer that is less than the number of 
equations. We will discuss this more in Section \4. 21 

The convergent path along the attractor is controlled by 
the parameters in the governing equations. When a control 
parameter i s slowly increa sed, the system exhibits a series of 
behaviours (|Scargle!ll992l ). This series is known as the 'route 
to chaos'. The behaviours usually unfold as: constant => ■ pe- 
riodic =>■ period two =>■ . . . chaos (Olsen & Dcen 1985|), as 
demonstrated in Fig. [3] where period two is a repeating path 
that travels twice around the attractor. As the parameter is 
increased, this behaviour continues to where the path cycles 
three, four, or more times to return to the same location, 
until it suddenly becomes chaotic, where the path will never 
repeat. 



3 LINEAR ANALYSIS 

We wish to search for chaotic and non- linear behaviour i n 
the spin-down rate presented in fig. 2 of iLvne et al.l (|20 lOh . 
There they isolated a subset of 17 pulsars with prominent 
variations in their frequency derivatives. These time series 
are ideal for non-linear studies because they directly relate 
to the dynamics of a pulsar. Before non-linear analysis can 
be done, we need to compensate for some their limitations. 



3.1 Mind the gap 

When dealing with large time-scales, such as those encoun- 
tered in astronomy, it is not always feasible to record data 
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Figure 2. Three dimensional visualization of two chaotic attractors. (a) The Lorenz attractor from t = to 100 with a = 10, p = 28, 
= 8/3 and x a = 0. y a = 10, z a = 10.2. (b) The Rossler attractor from i = to 500 with a = 0.2, b = 0.2, c = 5.7 and x a = -1.887, 
y = -3.5, z = 0.09789 
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Figure 3. The route to chaos for the Rossler equations 
a = b = 0.1 
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Figure 5. (a) Sine function: an example of a continuous function. 
(b) Random time series with the same number of data points as 
the sine wave. 



at regular intervals. This produces a times series that spo- 
radically samples a continuous phenomenon. When the time 
spacing between two points in the series is relatively small, 
little information is lost about the continuum inside that re- 
gion. If the spacing is large, more information is lost, which 
can cause a change in the structure of the data. 

To avoid such changes, we would like to analyze the 
largest section in the time series that best samples the phe- 
nomenon. We start by finding the statistical mode of the 
spacing, which gives us the step size that is the closest to 
being evenly sampled. We then compare this with the spac- 
ings in the series. If a gap is greater than three times the 
mode spacing, we assume that this region has significant in- 
formation loss. The time series is now broken into several 
sections that are separated by these large gaps, an example 

We extract the longest 



of which can be seen in Fig. 4(a 



section of data which is then normalized on both axes for 



more efficient computing, as seen in Fig. 4(b) 



3.2 Turning point analysis 

We want to ensure that the new time series is depicting 
a phenomenon and is not solely a consequence of noise. If 
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Figure 4. (a) The time series recorded for PSR B0740— 28. The vertical red bars highlight the gaps that are greater than three times 
the mode spacing, (b) The largest section for PSR B0740— 28 after being normalized on both axes. 
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Figure 6. Illustrations showing the combinations of data points 
for which u\ < u 2 < 113. See text for details. 



the time series samples a continuous function well, then we 
should expect a smooth curve with a large number of data 
points between maxima and minima, or what are known as 
turning points. If the time series samples a random distribu- 
tion, then we should expect several rapid fluctuations which 
will cause a large number of turning points. We can see an 
example of this in Fig. [5] where the well sampled continuous 
sine wave has very few turning points, while the random se- 
ries, with the same number of data points, has considerably 
more. 

In order to compare a series to noise, we need to have 
some expectation of the number of turning points that 
should be in a random series. It takes three consecutive data 
points to create a turning point. When randomly sampling, 
it is virtually impossible to have two values that are ex- 
actly the same. Therefore, the values of th e three data points 
will h ave the relationship, u\ < u 2 < 113 l|Kendall fc Stuartl 
1966). These values can be rearranged in six different combi- 
nations, as seen i n Fig.^l were we find fou r that will produce 
a turning point jKendall fc StuartJll966ft . This means that 
any enclosed data point will have a 2/3 probability of being 
a turning point. The first and last data points in a series 
of size n are the only unenclosed points. Thus, there are 
(n — 2) possible locations where a turning point can occur. 
This leads to an expectation value 



2 



90 



(4) 



Mr = 2) 



(3) 



If the number of turning points is greater than the ex- 
pected value, the series is fluctua ting more rapidly than ex - 
pected for a random time series jBrockwell fc Davis! Il99ri ). 
On the other hand, a value less than the expected value in- 
dicates an increase in the number of data points between 
each turning point. 

To be confident that the signal is not a result of random 
fluctuations, we will only keep the series that have a total 
number of turning points that are five standard deviations 
less than the expectation value of a series of the same size. 
If a time series has a total number of turning points greater 
than this value, it is indistinguishable from noise and will 
not be used. This does not necessarily imply that there is 
not a phenomenon present, rather that the data may not 
sufficiently sample the phenomenon. 

In Table [1] we have listed the values for the time series 
that was extracted from the corresponding pulsar. Unfor- 
tunately, the last four pulsars failed to pass the detection 
threshold, and are excluded from any further analysis. 

3.3 Cubic spline 

Most signal processing techniques assume that a time series 
is evenly sampled, and when the series is spaced randomly 
these algorithms severely increase in complexity. Therefore 
we would like to form an evenly sampled series that has 
the same structure as our own. Since we have selected a 
section that is nearly evenly spaced and we are reasonably 
confident that a signal is resolved, we can use a cubic spline 
interpolation to approximate the structure in between data 
points. After the data have been splined, we then resample 
with a step size equal to the statistical mode of the original 
spacing. This is done with the intent of limiting a significant 
increase in the time resolution. 

3.4 Fourier transform 

Having more data points will be beneficial in upcoming anal- 
yses, but caution must be used to avoid introducing struc- 
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Pulsars 


T 


n 


HT 




# of cry's away 
T is from fij< 


B1828-11 


56 


259 


171.33 


6.76 


-17.06 


B0740-28 


105 


270 


178.67 


6.90 


-10.67 


B1826-17 


34 


123 


80.67 


4.64 


-10.05 


B1642-03 


53 


151 


101.33 


5.20 


-9.29 


B1540-06 


19 


74 


48.00 


3.58 


-8.10 


J2148+63 


25 


81 


52.67 


3.75 


-7.37 


B0919+06 


31 


90 


58.67 


3.96 


-6.99 


B1714-34 


7 


39 


24.67 


2.57 


-6.87 


B1818-04 


23 


73 


47.33 


3.56 


-6.84 


J2043+2740 


12 


41 


28.00 


2.74 


-5.84 


B1903+07 


30 


79 


51.33 


3.70 


-5.76 


B0950+08 


29 


76 


49.33 


3.63 


-5.60 


B1907+00 


32 


80 


52.00 


3.73 


-5.36 


B2035+36 


14 


43 


27.33 


2.71 


-4.93 


B1929+20 


11 


33 


20.67 


2.35 


-4.11 


B1822-09 


56 


108 


70.67 


4.34 


-3.38 


B1839+09 


13 


31 


19.33 


2.28 


-2.78 



Table 1. The 17 pulsars listed in order of their distance away from 
fiT. Here T is the total number of turning points in the extracted 
series, n is the total number of data points in that series. (j,t is 
the expected number of turning point for a random series of size 
n, and ctt is the standard deviation of the expected value. The 
dashed line marks the five standard deviation threshold. 

tures or signals that are not present in the series. If we take 
the Fourier transform of the data, this will give us an am- 
plitude and a phase at each frequency below the Nyquist 
frequency. It is important to note that we receive both am- 
plitude and phase only when the data are evenly sampled. 
These can in return be used to construct an approximate 
time function for the series. This function can then be sam- 
pled at any rate without adding any significant non-pre- 
existing frequencies as described below. 

3.4-1 Noise reduction 

Even though the turning point analysis convinced us that 
our series is not governed by noise, this does not mean that 
there is no noise in the data. Before generating the new time 
series, we take the opportunity at this point to perform noise 
reduction. 

The simplest method is a low pass filter, but a decision 
needs to be made as to where to place the cut-off frequency. 
One can view our time series as the addition of two series: 
a signal series plus a noise series of equal size. From our 
turning point analysis, we know the expected number of 
total turning points for the noise series, and on average we 
should expect the minima and maxima to be equally spaced. 
Therefore, on average the noise should contain a frequency 



were At to tai is the total time recorded. We then extend the 
frequency to a value corresponding to two more standard 
deviations below the expected number. Our data have been 
normalized in time, so that they only cover a single time 
unit. This leads us to set an angular cut-off frequency 

CJmax = Tr(fi T ~ 2a T )- (6) 
We are now able to approximate a time function by 



doing a reverse Fourier transform up to the cut-off frequency. 
Since the Fourier series is a sum of periodic functions, the 
new time series will, by definition, repeat back upon itself. 
This can cause noticeable errors towards the beginning and 
end of the time series, as seen in Fig. [7] and to avoid this we 
will not sample within the first and last five percent of the 
time recorded. We then generate the time series by evenly 
evaluating the time function until we have 5,550 data points. 
This leaves us with approximately 5,000 data points within 
the desired range. 



4 NON-LINEAR ANALYSIS 

For the non- linear analysis we use a combination of the 
TISEAN software package presented bv lHegger et alj (|l998T ) 
and our programs written in the matlab programing lan- 
guage. We intend to have the matlab programs available to 
the reader on MathWorks file exchange website to help aid 
in the understanding of the algorithms. 

4.1 Attractor reconstruction 

With only one observable, it seems that we are unable to 
reproduce an attractor, but surprisingly the dynamic series 
contains all the information needed for its reconstruction. 

4-1.1 Method of delays 

We can u se embedding theor ems developed bv lTakensl (| 1 98 lh 
and by ISauer et all i|l99ll ) to reconstruct the attractor. 
These theorems state that a series of scalar measurements 
of a dynamical system can provide a one-to-one image of a 
vecto r set, the strange at tractor, through time delay embed- 
ding (Hcg ger~et al.lfl998T ). Each element in the vector is the 
scalar measurement S(t) at a different time as follows 

x(i) = [S(t), S(t + t), S{t + 2t), . . . , S(t + (m - 1)t)}. (7) 

The number m of ele ments in the vector is said to be the 
embedding dimension (|Hegger et al.|[l99^ ') , while r is a time 
delay. 

4-1.2 Time delay 

It soon becomes evident that picking the proper time delay 
is crucial. If the time delay is too small, then each element 
will be very close in value, forming a tight cluster. If the 
delay is too large, then the elements are unrelated and the 
attractor information is lost. If we were simulating a solu- 
tion to the non-linear equations, we would need to find a 
time-scale where linear effects are dominant and to make 
sure that our step size was within this range. We would like 
to do the same thing but for the time delay, because the 
orthogonal axes will differ by a dynamically linear trigono- 
metric function. 

There are a wide range of algorithms that are used to 
find this appropriate time delay, but the simplest on es fail 
to account for non- linear effects (|Hegger et al.l B-998). The 
algorithms that do account for these effects are not very 
intuitive. We decided to use ad aptations of turbulent fl ow 
techniques which can be seen in iMathieu fc Scottl |2000l ). 
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Figure 7. (a) The Fourier transform amplitudes of PSR B1642-03. The vertical red line marks the cut-off frequency, (b) Top: The pulsar 
data (red) overlaid with Fourier time function with the same cadence (blue). Bottom: The difference of the two curves. Note that at the 
end of the series the magnitude of the difference increases sharply. 





(a) 



(b) 



Figure 9. (a) The Lorenz attractor XY plane under same condition as Fig. |2(a)| (b) The reconstructed attractor from only the X values 
with t = 0.1345, which was estimated using the autocorrelation coefficients. 



One way of estimating this linear range is to use the au- 
tocorrelation coefficients. We are interested in the structure 
of the fluctuations and would like to have a zero average 
signal s(t). This is easily done by subtracting the time av- 
erage of the series from each scalar measurement. We then 
generate the autocorrelation coefficients for this new signal, 
defined as 



p(At) = 



< s(t)s(t + At) > 
< s(i) 2 > 



(8) 



This function will generate a curve that will start at one and 
then taper to zero, seen in Fig. [8] 

We can estimate the linear region by fitting a parabola 
to the small At region of p(At). The positive root of this 
parabola is our estimate, A, for a linear time-scale. This 
means that time-scales up to A should be dominated by lin- 
ear effects, but to make sure we are well into this region, 
we set the time delay to half of this value. We can see how 
well this estimate works in Fig. O where the topology of the 
attractor has been conserved. 

When these techniques are applied to the pulsar time 



series, a similar topology appears among the best-sampled 
pulsars, which is seen in Fig. 1101 This seems to suggest that 
these are different paths along the same attractor, demon- 
strating a route to chaos. One could see PSR B1540— 06 as 
periodic behaviour, B1828— llQas period two behaviour, and 
B1826-17 and B1642-03 as chaotic. To be sure that this 
is truly the case, we need to measure the dimension of each 
topology. 

4.2 Measuring dimensions 

We often see dimensions as the minimum number of coor- 
dinates needed to describe every point in a given geometry 
(Strogatz 1994). For example, a smooth curve is one dimen- 
sional because we can describe every point by one number, 
the distance along the curve to a reference point on the curve 



A linear trend was removed in order to make the time series 
stationary, because we are solely interested in the fluctuations. 
When B1828— 11 is referenced in the rest of the paper, it is implied 
that this linear trend has been removed. 
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Figure 8. An ideal example of autocorrelation coefficients. The 
blue curve is the autocorrelation coefficients. The red dashed 
curve is a parabolic fit to the autocorrelation on small time-scales. 
A is the time-scale estimate where linear relationships are believed 
to be dominant. 
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Figure 10. The topology of the pulsar time series after being 
embedded in three dimensional space. 



ilStrogatzlll994h . This definition fails us when we try to ap- 
ply it to fractals, and we can see why by examining the von 
Koch curve. 

Shown in Fig. 1111 the von Koch curve is an infinite limit 
to an iterative process. We start with a line segment 5b and 
break the segment into three equal parts. The middle section 
is swung 60 degrees and another section of equal length is 
added to close the gap to form segment Si. This process is 
repeated to each line segment to produce the next iteration. 

With each iteration, the curve length is increased by a 
factor of |. Therefore, the total length at an iteration n will 



Figure 11. Iterations leading to the von Koch curve. S n denotes 
the segment on iteration n and the far right column is the total 
length of the corresponding segment. The von Koch curve is when 
n — > oo. 



be L n — Lo x (§)"• When this is iterated an infinite amount 
of times, it produces a curve with an infinite length. Not 
only that, there would be an infinite distance between any 
two points on the curve. 

We would not be able to describe every point by an 
arc length, but we would be able to describe them with two 
values, say the Cartesian coordinates. Our original definition 
would define this as two-dimensional. However, this does not 
make sense intuitively, since this is not an area. Therefore 
this is something in between, one plus some fraction of a 
dimension. This is known as a fractal dimension. 



4-2.1 Correlation dimension 

There are several different algorithms that are used to mea- 
sure this fractal dimension, but nearly all depend on a 
power-law relationship. The most widely used is the correla- 
tion d imension, first introduced bv lGrassberger fc Procaccial 
(1983ft . 

It is calculated by creating a test sphere of radius R 
centred on a data point located at x. The number of data 
points inside the sphere is then counted, N X (R). The radius 
is slowly increased. As it increases, t he number of da ta points 
in the sphere grows as a power-law (Strogat j|l994f ) 



N X (R) oc R d 
d(lnN x ) 



d = 



d(ln R) 



(9) 



Here d is referred to as the pointwise dimension at x. Fig. 
[12] shows this relationship with familiar geometries. 

Due to fluctuations in the sampling density, the point- 
wise dimension can vary depending on where x is located. 
In order to produce a more self consistent measurement, an 
average of N x (R) is taken over all data points. This average 
is known as the correlation sum, C(R), and is often written 



C(R) 



N(N — 1) 



£ ^(fl-llxi-xj 



(10) 



i=l j=i-\-l 
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each of these situations. 



where H is the Heaviside function, H(x) — if x and 
H(x) — 1 if x > 0. N is the total number of locations. || | 
is the magnitude of the vector such that ||v|| = (v • v) 1 / 2 . 

This correlation sum will have the same type of power 
relation as N X (R). The exponent for C(R) is then properly 
named the correlation dimension. This is defined as 



D 



lim lim . 
N-+oo r->o d(\nR) 



(11) 



We seldom have the luxury of infinite sample size with in- 
finitely small resolution, and therefore different behaviours 
occur in the correlation sum. 

To estimate D, one would plot InC versus \nR. If the 
power relationship held true for all R we would see a straight 
line with a slope of D. However, due to the finite size of 
the attractor, all points could lie within a large enough R. 
We can see in Equation I1UI that this would cause C(R) to 
converge towards a value of one. At low enough R we will 
reach a resolution limit, where only the centre data point 
is inside the test sphere. This causes the C(R) to converge 
to zero. Therefore the power-law will only hold true over an 
intermediate scaling region. 

We could safely avoid this resolution limit if 10 points 
were in each of our test spheres. This would correspond to 
C(R) = tS, which we use as our lower limit of the scaling 
region. To set an upper limit, we say that if 10% of all points 
are in the each test sphere, we would start to approach the 
scale of the attractor. Therefore, we set our upper bound at 
C(R) = 0.10 to close off a rough scaling region. 

With this definition of dimension, one can see how there 
is a transition between integers. We can increase the num- 
ber of data points of a line within R by simply folding that 
line inside a test radius, perhaps like Si in Fig. 1111 In or- 
der to have a power-law relationship across all sizes of R, 
we need to have a similar fold on all scales, known as self 
similarity. If we were to steepen the angles of each fold, this 
would increase the number contained and also increase the 
dimension. We can do this until the folds are directly on top 
of one another to 'colour in' a two dimensional surface. We 
could then repeat the process by folding the two dimensional 
surfaces to transition to three dimensions. Folding like this 
is seen as the underlying reason why strange attracto rs have 
their fract al dimensions and is explo r ed in depth by ISmalel 
l| 19671 ) and lGrassberger fc Procaccial (jl983r i. 



Figure 13. Space time separation plot, generated using the 
TISEAN package, for PSR B1642— 03 with m = 3. The contours 
indicate the percentage of locations within a distance at a given 
time separation. The different curves correspond to 10% increases 
starting with the lowest curve at 10% to the top curve at 90%. 



4-2.2 Theiler window 

If the attractor was randomly sampled we would be ready 
to measure the dimension, but the solutions are one dimen- 
sional paths on a multi-dimensional surface. Therefore, we 
have two competing dimension values which will corrupt our 
measurement. What we wish to do is to isolate only the at- 
tractor dimension, the geometric correlations, and remove 
the path relationships, the te mporal correlatio ns. 

As difficult as this sounds, iTheilerl |l986J) introduced a 
rather simple remedy. We exclude the points within a time 
window, wx At a , around the centre of the test radius. This 
is known as the Theiler window and changes the correlation 
sum in the following way 



C(R) = 



(N-w)(N-w-l) 



J2 E m-\\*i- 



i= 1 j=i-\-w-\-l 



(12) 

We can see that picking the appropriate w is nontrivial. 
If we pick a window that is too small we fail to remove 
the temporal correlations, and if the window is too large 
it would significantly reduce the accuracy of the geometric 
measurement. 



4-2.3 Space time separations 

In order to estimate a safe value for the Theiler window, 
IProvenzale et all (Il992l ) introduced the space time separa- 
tion plot. It shows the relationship between the spatial and 
temporal separations, by forming a contour map of the per- 
centage of locations within a distance, for a given time sep- 
aration. 

An example of this is seen in Fig. 1131 here we used the 
stp program in the tisean package. We can see that the lo- 
cations with small time separations are close to one another. 

2 Here At a is the step size of the time series, and w is an integer 
value. 
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As the time separations increase, so do their spatial sepa- 
rations, until they reach some asymptotic behaviour. The 
substantial fluctuations at larger times are contributed t o 
the cycle period of the attractor (|Kantz Sz Schreiberll20o4 ). 
Tempora l correlations are present until the contour curves 
saturate ijKantz fc Schreiberl2004h . In Fig.[T3]this transition 
seems to occur around a time separation of 40 steps. This is 
a similar quantity that appears in all of our data sets, but 
to be safe we chose a more conservative value of w = 100 
steps for the rest of our analysis. 

4-2.4 Embedding to higher dimensions 

We are now ready to measure the dimensions of the topolo- 
gies that were created in Section 14.11 If these topologies 
are true geometric shapes, their dimensions will not change 
when placed in a higher dimensional space. This means that 
as long as our embedding dimension, m, is greater than 
the attractor dimension, our correlation dimension measure- 
ment will be constant. On the other hand, if we were to em- 
bed a random distribution, it would be able to occupy the 
entire space and would have a correlation dimension equal 
to the embedding dimension. 

Because of this behaviour, we can embed our time series 
to higher and higher dimensions to see if it will plateau to a 
constant. We cannot do this forever because the number of 
location vectors that could be formed drops off as 

N = n s -T*(m-l), (13) 

where n s is the total number of scalar values in the original 
time series and r* is the number of time steps, At s , corre- 
sponding with the time delay, r. This drop off causes a slight 
reduction in our statistical accuracy as we continue to higher 
dimensions. Therefore, we limit our embedding dimensions 
to between 2 and 10. 

Though we have picked a rough scaling region in Sec- 
tion l4.2.Tl C ~ — to 0.10, this does not mean that the true 
scaling region extends through this entire range. Therefore, 
we sweep the region with a window size from a third of the 
total logarithmic range to the entire range, searching for the 
minimum deviation in the correlation dimension measure- 
ments over the embedding dimensions. This ensures that we 
are picking the 'flattest' possible section over a significant 
scale. 

When this technique is applied to the four pulsars from 
Fig. 1101 two pulsars plateau very nicely. PSR B1828— 11, 
seen in Fig. 1141 averages to a correlation dimension of 
2.O6±O.O30, and PSR B1540-06 averages to 2.50±0.09, for 
embedding dimensions greater than 3. The other two pul- 
sars fail to converge to a constant. We attribute this non- 
convergence to the sparse coverage, where each pass around 
the attractor is too far apart to get a proper measurement. 

4-2.5 Surrogate data 

We want to guarantee that any plateaus are due to geometric 
correlations and cannot be produced by random processes. 
Therefore, we would like to test several different time series 

3 The error calculation is presented in appendix [51 




Figure 14. (a) InC versus laR for PSR B1828-11. The dashed 
horizontal blue lines mark the rough scaling region, C = — to 
0.10. The solid horizontal red lines mark the flattest scaling region 
that is larger than a third of the rough scaling region, (b) The 
flat blue line is the correlation dimension measurements, within 
the scaling in (a), versus embedding dimension. The line at 45° is 
what a purely random time series would be in that given m. The 
sloped red line is the mean and standard deviation of 10 surrogate 
data sets (see text). 

with the similar mean, variance, and autocorrelation func- 
tion as the original data set. These data sets are known as 
surrogate data. 

The idea of s urrogate data was first introduce in 
iTheiler et af] d 19921). but we use an improved version that 
was presented in Ischreiber fc Schmit3 (|l9991 1 for our sur- 
rogate data sets. They start by shuffling the order of the 
original time series, and then take the Fourier transform of 
this shuffled series. Keeping the phase angle of the shuf- 
fled set, they replace the amplitudes with the Fourier am- 
plitudes of the original times series. Then they reverse the 
Fourier transform, and crea te the desired surrogate data set. 
ISchreiber fc Schmitj (|l9991 ) iteratively do this until changes 
in the Fourier spectrum are reduced. Fortunately, this is 
all done in the pro gram surrogates in the tisean package 
l|Hegger et al.| [l998). which we used for our analysis. 

We form 10 surrogate data sets for each pulsar and run 
them through the same non-linear analysis as the post pro- 
cessed time series. We then find the average and standard 
deviation for the surrogates' correlation dimension for each 
embedding dimension to trace out the region were we should 
expect other surrogates to lie. We are confident that the 
original time series is not due to a random process if its cor- 
relation dimension measurements lie away from this region. 

As seen in Fig. [14] PSR B1828-11 is well outside the 
surrogate region, suggesting that it is a true dimension mea- 
surement. The measurements for PSR B1540— 06, seen in 
Fig. 1151 marginally misses this boundary, but because this 
pulsar is highly periodic, the Fourier spectrum is dominated 
by a single frequency. This restricts the complexity of the 
surrogate to where little change in dimension is expected. 
We will see more on this in Section T4. 2. 61 

4-2.6 Correlation Benchmark Testing 

Though each part of the algorithm works independently, we 
want to ensure that it all works together correctly. Therefore 
we need to run through the process with a time series from 
known attractors under similar conditions to our original 
data, in order to see if we receive similar dimensions. For this 
analysis we chose the Lorenz and Rossler attractors, whose 
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(a) (b) 

Figure 15. (a) The amplitudes of the Fourier series versus nor- 
malized frequency for PSR B1540— 06. The vertical red line marks 
where we set the cut of frequency, (b) The blue error bars show 
the correlation dimension measurement, and the red error bars 
show the surrogate data region for PSR B1540— 06. 



dimensions have been do cumented as being 2.07 ± 0.09 and 
1.99 ± 0.07 respectively (|Sprottll200 $ ). 

We first generate the x component time series for a cor- 
responding attractor, making sure that our initial conditions 
are well on the attractor. Then, we extract a section from 
this time series out to the same number of turning points 
as used in the pulsar time series. This extracted section was 
then normalized on both axes, as we did in Section \3. II We 
then resample the normalized time series by applying a cu- 
bic spline to the normalized times for the pulsar time series. 
This ensures that the new time series has the same time 
spacing and number of turning points as the pulsar time se- 
ries. This new time series is then run through the algorithm 
to see if the proper dimension is achieved. 

An example from the correlation dimension measure- 
ments are seen in Fig. 1161 where we can see that due to the 
simplicity of the Rossler attractor in the frequency domain 
the surrogate data region only deviates by a small amount. 
We refer to this as a borderline detection. We also classify de- 
tections as being inside or outside based on their proximity 
to the surrogate region. Regardless, the correlation dimen- 
sions for both attractors plateau for all of the pulsar time 
series. 

The results from all of the series are listed in Table [2] 
We can see that the algorithm does rather well considering 
we start with less than 300 data points. Although the algo- 
rithm is not refined enough to pinpoint the dimension, it is 
sufficient to determine the number of governing variables. If 
we round the dimension measurement to the nearest inte- 
ger and then add one we will receive the proper number of 
variables, for any series that is either outside or borderline 
to the surrogate region. 



4.3 Measuring the butterfly effect 

Though our correlation dimension measurement can narrow 
the number of governing variables, it is not definitive enough 
on its own to say that our attractors are chaotic. Therefore, 
we need to search for other signs of chaos in our topologies 
to ensure that they are indeed strange attractors. 




Figure 16. Benchmark results for PSR B1828— 11. Left column 
is the Lorenz attractor results. Right column is the Rossler attrac- 
tor results. Top row are the reconstructed attractors. Middle row 
are the Fourier transforms where the vertical red line marks the 
cut off frequency from Section l3.4.1l Bottom row are the correla- 
tion dimension measurements (blue error bars) with the surrogate 
region (red error bars). 



Lorenz Rossler 





Correlation 


Surrogate 


Correlation 


Surrogate 


Pulsar 


Dimension 


Region 


Dimension 


Region 


B1828-11 


2.20 ±0.06 


O 


1.96 ±0.02 


B 


B0740-28 


2.11 ± 0.03 


O 


2.16 ±0.08 


B 


B1826-17 


2.00 ±0.10 





1.96 ±0.04 


B 


B1642-03 


2.05 ± 0.07 





2.00 ±0.10 


B 


B1540-06 


1.90 ±0.20 





2.30 ±0.10 


B 


J2148+63 


2.33 ± 0.05 





2.13 ±0.02 


B 


B0919+06 


1.90 ±0.10 





2.00 ±0.10 


1 


B1714-34 


1.85 ± 0.01 


1 


2.87 ±0.07 


I 


B1818-04 


1.70 ±0.05 





2.26 ±0.08 


B 


B2044+2740 


1.83 ± 0.03 





2.10 ±0.10 


O 


B1903+07 


2.11 ± 0.09 





2.00 ±0.20 


I 


B0950+08 


1.68 ± 0.09 





2.00 ±0.20 


B 


Average 


2.00 ±0.20 




2.00 ±0.10 





Table 2. The weighted average for the correlation dimension for 
m > 3 for the corresponding pulsar and attractor. The surrogate 
region column marks weather the plateau line was either inside 
(I), outside (O), or borderline to (B) the surrogate data region. 
The average row is the weighted average of the correlation dimen- 
sion measurements that were either marked as O or B. 
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4-3.1 Lyapunov exponent 

As mentioned in Section[2l the butterfly effect is only present 
in chaotic systems. We should only see an exponential di- 
vergence in nearby locations over time if the system is 
non-linear. The exponent of this increase is a character- 
istic of the system and quantifies the strength of chaos 
(Kan tz fc Schreibeij[2003 ). This is known as the Lyapunov 
exponent. 

There are as many Lyapunov exponents as there are 
axes. Here we are only interested in measuring the max- 
imum exponent because it gives us the most information 
about our system. We start by looking at the distance be- 
tween two locations, 5q = ||xt x — x t 2 ||j then recording this 
separation over time, <5a* = Hxtj+At — xt 2 +At||- The Lya- 
punov exponent, A, would be 



(14) 



Since the separation between two points cannot be greater 
than the size of the attractor itself, Equation [TJ will only 
hold true for At values where 5 At is smaller than the attrac- 
tor. 

If A is positive, the system will be highly sensitive to 
initial conditions and therefore chaotic. If A is negative, the 
system would eventually converge to a single fixed point. If 
A is zero then this would represent a limit cycle where the 
path keeps repeating itsel f and is said to be marginally stable 
|Kantz fc Schreiberil2004h . 



In actuality, the separations do not grow everywhere 
on the attractor and locally they can even shrink, and with 
contributions from experimental noise it is more robust to 
use the average to obtain the Lyanpunov exponent. 

The algorit hm for this averaging was introduced in - 
dependently by iRosenstein et al.l j 19931 ) and iKantj (|l994h . 
One first picks a centre point located at xt and takes note 
of the data points within a radius of e. This is known as 
the neighbourhood. For a fixed At, the mean $At is cal- 
culated across the whole neighbourhood. The logarithm of 
this mean distance is then averaged over all points in the 
attractor. Therefore, one needs to compute 



1 tN 

S(At) = - In 



\U(xt )\ 



E 



Xt + At — Xt+At 



*tet/(x to ) 



(15) 

where C/(xt ) is the neighbourhood centred on xt . For our 
analysis we chose an e that is twice as large as our aver- 
age distance to the nearest neighbour, while accommodat- 
ing for the Theiler window. If the exponential relationship 
is present, we should see an overall linear behaviour when 
S(At) is plotted with respect to At. In order to conserve 
statistical accuracy we do not look beyond time-scales that 
are half the total time. 

Similar to the correlation dimension, the Lyanpunov ex- 
ponent is invariant to the number of embedding dimensions 
as long as m > d. Therefore we calculate S(At) from m = 2 
to 10 to ensure that the slopes in the linear regions remain 
constant. When noise is present in deterministic systems, 
it causes a process similar to diffusion where t?At expands 
prop ortionate to \f~A~i on small scales (|Kantz fc Schreiber] 
2(1111:. This causes S{At) to have a ^ln(At) behaviour, 
which produces a steep increase over short time intervals. 

When this is applied to PSR B1828-11, seen in Fig.[T7l 



three distinct regions appear. On small scales, we can see a 
convergence re gion due to a combi nation of a noise floor and 
non- normality (Hcg ger et al.lfl998l ). Because of these effects, 
it takes the neighbourhood a while to align in the direction 
of the largest Lyapunov exponent. Once the neighbourhood 
has converged, we see a similar linear behaviour across all 
the embedding dimensions. This continues until the curves 
starts to saturate to a value on the scale of the attractor for 
that particular embedded dimension. 

The positive slope of the linear region suggests that 
B1828— 11 is a chaotic system with a maximum Lyapunov 
exponent of (4.0 ± 0.3) x 10 -4 inverse days. When the same 
procedure is done to a surrogate data set of B1828— 11, no 
linear region is observed. Once the surrogate neighbourhood 
passes the convergent region it directly saturates around 
a constant value. This constant saturation would suggest 
a Lyapunov exponent of zero, which is what is expected 
for a linear dynamical equation. The surrogate's strikingly 
different behaviour would imply that the linear region in 
B1828— 11 is a real detection and is not a consequence of a 
random process. 

When this is applied to the other pulsars there are no 
other definitive detections. Either the convergent region is 
dominant causing a direct saturation, or the total time of 
the measurement is too small to positively state a separation 
between the convergent and linear regions. We will expand 
on these ideas more in the following section. 



4-3.2 Lyapunov Benchmark Testing 

Again we want to ensure that the algorithm is working in its 
entirety to produce the known values of the Lyapunov ex- 
ponent. Therefore we follow the same procedure as Section 
14.2.61 matching the same number of turning points and sam- 
pling spacing of the pulsar. Because PSR B1828— 11 is our 
only detection, we concentrate our attention on its bench- 
mark results. 

The Lyapunov exponent changes with different param- 
eters in the governing equations. Because of this we use the 
most common chaotic parameters for the Lorenz and Rossler 
attractorfl Under these conditions, the maximum Lya- 
punov exponent for the Lorenz attra ctor is A max — 0.9056 
and for the Rossler is A max ~ 0.0714 (|Sprottll200l '). 

The results for the benchmark testing of B1828— 11, 
seen in Fig. 1181 seems to produce conflicting results. The 
algorithm estimates a maximum exponent for the Rossler 
attractor to be a reasonable A max = 0.080 ± 0.003, but for 
the Lorenz attractor it estimates an unreasonable value of 
A max = 0.28 ± 0.05. 

The cause of this discrepancy can be seen in the sur- 
rogate results. There we can see that the convergence time- 
scales are affecting our estimates. For the Rossler attrac- 
tor, the time for the neighbourhood to converge is about 
5 time units, while the exponential behaviour will last on 
time scales of — - — ~ 14 units. Therefore, the Rossler at- 
tractor is outlasting the convergence time and will be able to 
demonstrate its exponential behaviour. The Lorenz attrac- 
tor with these parameter values is more sensitive to initial 
conditions and will only exhibit its exponential behaviour 



4 Lorenz: a = 10, r = 28, b = 8/3; Rossler : a = b = 0.2, c = 5.7 
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At At 

(a) (b) 

Figure 17. (a) S(At) versus time elapsed for PSR B1828— 11. Each curve is for a different embedding dimension from 2-10 starting 
at the lowest curve. The vertical red lines mark the region where least square regression was applied. This corresponds to an average 
maximum Lyapunov exponent of A max = (4.0 ± 0.3) X 10 — 4 inverse days, (b) S(At) for a surrogate data set for B1828— 11. 




Figure 18. Benchmark results for PSR B1828— 11. Left column is 
the Lorenz attractor results. Right column is the Rossler attrac- 
tor results. Top row are the S(At) for reconstructed attractors. 
Bottom row are the S(At) surrogate data for the appropriate 
attractors. 



on a time-scale of j-= — w 1 unit, which is on the same scale 
of its convergence time. Therefore the S(At) for the Lorenz 
attractor is not portraying its true behaviour. 

This convergence time is partly inherent to the system 
and partly due to noise. Because of this there are two ways to 
correct this behaviour. The first is to simply reduce the noise 
in the time series. The other way is to increase the sampling 
period of the measurements in order to start with smaller 
neighbourhoods which will give the exponential behaviour 
several orders of magnitude to appear. Unfortunately, we 
are unable to comfortably remove any more noise in the 
algorithm. Also, to reduce the neighbourhoods by orders of 
magnitude, we have to increase the number of data points by 



orders of magnitudes, which is currently beyond our system's 
capabilities. 

Because of this convergence time, our algorithm is sen- 
sitive to lower levels of chaos with smaller Lyapunov expo- 
nents and is not able to distinguish between higher levels of 
chaos and the convergence region. 



5 CONCLUSIONS 

By using a careful combination of turning point analysis, 
cubic splining, and Fourier transforms, we have constructed 
an algorithm that re-samples an unevenly spaced time se- 
ries without losing structural information. We have demon- 
strated this through an array of benchmark testing with 
known chaotic time series under similar conditions to a given 
pulsar time series. This testing has shown that there are no 
significant changes in the correlation dimension or the max- 
imum Lyapunov exponents, when it was detectable. 

Thes e techniques were applied to the pulsar spin-down 
rates from lLvne et al.l (|2010l l. where PSR B1828-11 exhibits 
clear chaotic behaviour. We have shown that the measure- 
ments of its correlation dimension and maximum Lyapunov 
exponents are largely invariant across embedding dimen- 
sions. This, combined with its strikingly different reactions 
compared to its surrogate data sets, has shown that the 
chaotic characteristics in this pulsar are not caused by ran- 
dom processes. 

The positive measurement of A max = (4.0 ± 0.3) x 
inverse days confirms that B1828— 11 is chaotic in nature. 
For a system of equations to be chaotic there needs to be 
a minimum of three dynamical equations with three gov- 
erning variables. The correlation dimension measurement of 
B1828-11, D = 2.06 ± 0.03, implies that there are a to- 
tal of three governing variables, meeting the minimum re- 
quirements for chaos. One governing variable is clearly the 
spin-down rate of the pulsar. At this time we can only spec- 
ulate on the other two. We know that the magnetic fields 
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of a pulsar can change its dynamics and change the pulse 
profile. It has also been shown that changes to the super- 
fluid interior can affect the long term dynamics of a pulsar 
|Ho fc Anderssonll2012l '). Therefore, these seem to be great 
candidates. Beyond this, the subject is still a mystery and 
needs to be explored with non-linear simulations. Regard- 
less of the model chosen, it is possible to perform some of 
the methods presented here, on the simulated data, to see if 
similar chaotic behaviours are present. 

Knowing that B1828— 11 is governed by three variables, 
the reconstructed attractor in Fig. [TO] would be an accurate 
depiction of its strange attractor. Because this is visually 
similar to the other attractors in Fig. 1101 we would find it 
peculiar if their dynamics were not somehow related. Unfor- 
tunately, the techniques in this paper were unable to confirm 
this relationship. If these pulsars continue to be observed 
with an increase cadence, this would improve the correlation 
dimension and Lyapunov exponent measurements, perhaps 
to the point where these similarities could be quantified. 
With these measurements and a working model, estimates 
for the parameters can be given, which will give us further 
insight into the interior and/or exterior of these pulsars and 
how this relates to their dynamics. 
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APPENDIX A: TURNING POINT VARIANCE 
DERIVATION 



This derivation is presented in \Kendall & Stuar t (1966) but 
is rewritten here for the reader's convenience. 

The number of turing points p can be seen as a sum- 
mation 



(Al) 



where Xi is equal to one if Ui+i is a turning point, and zero 
if not. We have already shown that the expectation value of 
the sum 



E(j>) = Y / E{X i ) = \{n-2). 



(A2) 



In order to calculate the variance we need to find the expec- 
tation of the square of the number of turning points 



(A3) 



Now, substituting these values into lA5l we find 
E(p 2 ) = %{n - 2) + %(n - 3) + ^(n - 4) + Un - 4)(n - 5) 



E(p 2 



3 V ' 6 
40n 2 - 144n + 131 
90 



10 



9 



Finally, using the definition of the variance 



4 = E(p 2 ) - E(pf, 



(A10) 



(All) 



we insert our results from Equation I A2I and I A10I to conclude 
that 

2 40n 2 - 144n + 131 4(n - 2) 2 



DO 



!) 



2 16n - 29 



(A12) 



This can be expanded to 

E(p 2 ) = E | X ? + 2 E XiX i+ i + 2 XiX i+ 2 

K.n — 2 n — 3 n — 4 

+ Y, XiX i+k \,k^Q,l,2, (A4) 

(n-4)(n-5) J 

where the suffixes to the signs indicate the number of 
terms over which summation takes place. We then have 

E(p 2 ) = (n-2)E(X 2 )+2(n-3)E(X l X l+1 )+2(n^4)E(X l X t +2 
+ (n-4)(n-5)E(X l X l +k). (A5) 
Since X 2 = Xi, we have 



APPENDIX B: CORRELATION DIMENSION 
ERROR CALCULATIONS 

The correlation sum, d, is the average of the pointwise mea- 
surements around the attractor for a certain test radius, Ri. 
Because of this, we can use the standard deviation of a mean 
as the uncertainty of the correlation sum for that particular 
radius, 



& pointwise 



(Bl) 



ViV - w 

where N is the number of location vectors and w is the 
Theiler window. 

Next we use standard propagation of error techniques, 
to form the first order estimate of the uncertainty for the 
natural logarithm of the correlation sum, 

Q~Cj 



(B2) 



E(X 2 ) = l 



(A6) 



For k > 2, Xi and Xi+t are independent, for they have no 
values in common. Thus 



Wit h in the desired scaling region, we follow the outline 
in IPressI ((2007) for linear regression of least-squares to es- 
timate the correlation dimension and its uncertainty. IPressI 
l|2007l ) breaks this calculation into several summations. We 
have rewritten these sums for this application as follows: 



E(XiXi+k) = E(Xi)E[Xi+k) = 



9' 



(A7) 



We still need to evaluate E(XiX i+1 ) and E{XiX i+ 2). 
For the first term to contribute, there needs to be two 
consecutive turning points. In order to do this, we need 
a minimum of four data points, and define their values as 
mi < u-z < U3 < U4. This will generate 4! or 24 permutations, 
but one would find that only 10 permutations will contain 
two turning points. This leads to an expectation value of 



E(XiX i+1 ) - 23 - 12- 



(A8) 



Similarly, for E(XiXi+2) five data points with five dif- 
ferent values are needed. This produces 5! or 120 permuta- 
tions, but only 54 will produce a turning point in the second 
and fourth locations. This leads us to 



54 9 
E{X t X w ) = — = -. 



(A9) 



and 



N C 



'in B? 



N c 



S\ n R = 



In Ri 



i=N ff l»C, 



SlnC = E 



= N rfnCi ' 



N C 

E n . 



(In Ri 



In Riln d 



N C 

i=N "ItkOi 



(B3) 



(B4) 



(B5) 



(B6) 



(B7) 
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where iVo and Nc is the minimum and maximum index of 
the correlation sum measurements for a particular embed- 
ding dimension within the desired scaling region. With these 
summations, we are then able to calculate the slope. A com- 
mon denominator value, 



A = SS lnR 2 — (SinR.) 2 , (B8) 
is used to calculate the correlation dimension, 

, SSlnRlnC — SlnRSlnC 

dm = ^ , (B9) 



and variance, 



°in = f - ( B1 0) 



for each embedding dimension. 

To average the correlation dimension across the em- 
bedding dimensions, we use a weighted mean, where each 
measurement is weighted based on the inverse variance. For 
simplicity, we use normalized weighting coefficients, 

of 

Wm = —To "* _ 2 ■ (Bll) 

With this weighting, the mean correlation dimension, 

10 

< d >= w mdm, (B12) 

m=Z 

and the variance, 

10 

°<d> = w m (d m - < d >) 2 , (B13) 

m=3 

are calculated. 



